DALS016-高维数据推断1-背景介绍

前言

这一部分是《Data Analysis for the life sciences》的第6章线性高维数据推断(Inference For High Dimensional Data)的第1小节，这一部分的主要内容涉及高通量数据的介绍，统计学的一些注意事项，相应的R markdown文档可以参考作者的Github。

高通量技术

高通量技术已经将基础生物学和生物学从数据贫乏型学科尾成了数据密集型学科。一个典型的例子就是那些对基因表达有兴趣的研究领域。基因表达是这样一个过程，将生命蓝图DNA复制到RNA，而RNA是蛋白质合成的模板，是构成生命的基石。在20世纪90年代，对基因表达数据的分析相当于在一张纸上画了几个点，或从标准曲线上提取了几个数字。随着微阵列等高通量技术的发展，这几个数据就突然变成了筛选成千上万的数字。最近，RNA测序技术又进一步增加了数据的复杂性。生物学家们从使用他们的肉眼或简单的方法来分析结果到每个样本有数千（现在则是数百万个）个数据的分析。在这一篇中，我们佬双重点关注高通量数量的统计推断。具体来说就是，我们会使用统计检验来检测不同组之间的差异，以及使用有意义的方法来量化不确定性问题。我们还公介绍探索性数据分析方法，这些方法会结合高通量数据进行使用。在后面的章节中，我们会研究聚类，机器学习，因子分析和多层级建模(multi-level)。

由于现在有大量可用的公共数据集，我们会使用这些公共数据集来作为我们的学习案例。尽量使用的是公共数据，但是当你学习了这些方法后，在使用了高通量数据技术的其它领域也会派上用场。这些高通量技术包括微阵列，NGS，fRMI与蛋白质谱等。

数据包

我们这里使用的数据文件已经在作者的Github上分享了，因此使用devtools包就可以直接下载，如下所示：

library(devtools)
install_github("genomicsclass/GSE5859Subset")

3张表

我们用来做为案例学习的这些数据多数都是使用高通量技术生成的。这些技术会检测数千个特征值(feature)。这些特征值可以是基因，可以是基因组的单碱基位置，基因区域，以及图像的像素密度（芯片数据）。每个特定的检测产品被都会被相应地一组特定的功能所定义（我的理解就是，某公司的特定产品，例如检测miRNA的芯片，它的检测结果就是一些表达值，而这些表达值代表了哪些miRNA，就还需要另外一个miRNA名称的数据文件），例如特定基因表达的微阵列数据被一组其检测的基因所定义。一个特定的研究通常会使用一种产品对几个实验单元进行检验，例如几个人的样本。最常见的就是人，不过实验产单元也可以是其它的东西，例如肿瘤组织的一部分。按照实验术语，我们通常称实验单元为样本，这个样本与前面讲的随机样本不是一个东西，后者是统计学术语。

因此，一个高能被实验通常由3张表定义：一个是高通量数据的检测值，第二与第三表是关于列与行的信息。由于这些数据集通常是由一系列的实验单元，以及一些固定的特征构成，因此，高通量数据可以储存为$n \times m$的形式，其中$n$表示实验单元（样本），而$m$表示特征值（通常就是基因），现在 我们来看一个案例：

library(GSE5859Subset)
data(GSE5859Subset) ##this loads the three tables
dim(geneExpression)

结果如下所示：

> dim(geneExpression)
[1] 8793   24

这个数据是RNA表达谱数据，它检测了24个人（实验单元）的8793个基因。对于大多数的统计分析来说，我们还需要实验样本的一些信息。例如，在这个案例中，最初收集的数据是用于比较不同种族群体的基因表达情况。但是，我们已经创建了这个数据集的一个子集用于说明样本的信息情况，如下所示：

dim(sampleInfo)
head(sampleInfo)
sampleInfo$group

结果如下所示：

> dim(sampleInfo)
[1] 24  4
> head(sampleInfo)
    ethnicity       date         filename group
107       ASN 2005-06-23 GSM136508.CEL.gz     1
122       ASN 2005-06-27 GSM136530.CEL.gz     1
113       ASN 2005-06-27 GSM136517.CEL.gz     1
163       ASN 2005-10-28 GSM136576.CEL.gz     1
153       ASN 2005-10-07 GSM136566.CEL.gz     1
161       ASN 2005-10-07 GSM136574.CEL.gz     1
> sampleInfo$group
 [1] 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0

关于样本信息中的一列是文件名，它可以让我们把这个表的行与表达谱的列连接起来，如下所示：

match(sampleInfo$filename,colnames(geneExpression))

结果如下所示：

> match(sampleInfo$filename,colnames(geneExpression))
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

最后，我们还有一张描述特征值的表格，如下所示：

dim(geneAnnotation)
head(geneAnnotation)

结果如下所示：

> dim(geneAnnotation)
[1] 8793    4
> head(geneAnnotation)
     PROBEID  CHR     CHRLOC SYMBOL
1  1007_s_at chr6   30852327   DDR1
30   1053_at chr7  -73645832   RFC2
31    117_at chr1  161494036  HSPA6
32    121_at chr2 -113973574   PAX8
33 1255_g_at chr6   42123144 GUCA1A
34   1294_at chr3  -49842638   UBA7

从上面的信息知道，geneAnnotation这个文件里含有ID，这个表的行可以与表达谱的行连接起来，如下所示：

head(match(geneAnnotation$PROBEID,rownames(geneExpression)))

结果如下所示：

> head(match(geneAnnotation$PROBEID,rownames(geneExpression)))
[1] 1 2 3 4 5 6

实际应用

现在假设我们已经拥有了一批高通量数据，这个数据是检测了两组中一些人的基因表达情况。现在我们的问题是：在这两组人中，哪些基因的表达均值存在差异？如果是一个基因，计算起来很简单，我们手工计算就行，但是，高通常数量通常都是几千个基因，显然手工算不太现实，但是寻找差异基因的方法跟前面我们提到的统计推断原理一样。例如，我们可以使用t检验或其它检验来寻找差异基因。

随机变量的p值

这里我们要强调的一个概念就是，p值是一个随机变量。为了能说明这一点，我们来看一个案例。在这个案例中，我们会使用CLT近似的原来生成大量的样本，并通过t检验来计算其p值。然后我们的p值被定义为，正态分布的随机变量的绝对值大于观察到的t检验r值的概率，被称为$Z$（书中原话是：Then our p-value is defined as the probability that a normally distributed random variable is larger, in absolute value, than the observed t-test, call it Z.），因此，双尾检验的p值公式如下所示：

$$p = 2 \{ 1 - \Phi(\mid Z \mid)\}$$

在R中则为：

2*(1-pnorm(Z))

因为$Z$是一个随机变量，并且$\Phi$是一个确定的函数，$p$因此也是一个随机变量。我们现在使用Monte Carlo模拟来显示一下$p$值的变化，如下所示：

set.seed(1)

dir <- system.file(package = "dagdata")
filename <- file.path(dir,"extdata/femaleControlsPopulation.csv")
population <- read.csv(filename)
population <- unlist(population) # turn it into a numeric
N <- 12
B <- 10000
pvals <- replicate(B,{
  control = sample(population,N)
  treatment = sample(population,N)
  t.test(treatment,control)$p.val
})
hist(pvals)

从直方图中我们可以看到，p值属于均匀分布。事实上，我们可以从理论上证明，当零假设为真是，情况确实如此。当我们使用CLT时，我们的零假设$H_{0}$意味着我们的检验统计量Z服从均值为0，SD为1的正态分布，因此：

$$p_a = \mbox{Pr}(Z < a \mid H_0) = \Phi(a)$$

这也说明：

$$\begin{align*} \mbox{Pr}(p < p_a) &= \mbox{Pr}[ \Phi^{-1}(p) < \Phi^{-1}(p_a) ] \\ & = \mbox{Pr}(Z < a) = p_a \end{align*}$$

这其实就是均匀分布(uniform distribution)的定义。

千次检验

在上面的案例中，我们有两组数据，分别用0和1表示，如下所示：

library(GSE5859Subset)
data(GSE5859Subset)
g <- sampleInfo$group
g

结果如下所示：

> library(GSE5859Subset)
> data(GSE5859Subset)
> g <- sampleInfo$group
> g
 [1] 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0

如果我们对某个特定的基因感兴趣，例如第25行的那个基因，那么我们就可以简单地使用t检验就行，为了计算p值，我们需要先看一下这个基因数据是否服从t分布，因此我们可以使用qq图来绘制，如下所示：

e <- geneExpression[25,]
e
library(rafalib)
mypar(1,2)
qqnorm(e[g==1])
qqline(e[g==1])
qqnorm(e[g==0])
qqline(e[g==0])

结果如下所示：

> e
GSM136508.CEL.gz GSM136530.CEL.gz GSM136517.CEL.gz GSM136576.CEL.gz GSM136566.CEL.gz GSM136574.CEL.gz GSM136575.CEL.gz 
        10.34172         10.64662         10.69198         10.96155         10.24889         10.74414         10.22206 
GSM136569.CEL.gz GSM136568.CEL.gz GSM136559.CEL.gz GSM136565.CEL.gz GSM136573.CEL.gz GSM136523.CEL.gz GSM136509.CEL.gz 
        10.49094         10.57574         10.47758         10.43493         10.46440         10.58358         10.39607 
GSM136727.CEL.gz GSM136510.CEL.gz GSM136515.CEL.gz GSM136522.CEL.gz GSM136507.CEL.gz GSM136524.CEL.gz GSM136514.CEL.gz 
        10.75837         10.37383         10.57540         10.51203         10.39952         10.70539         10.77077 
GSM136563.CEL.gz GSM136564.CEL.gz GSM136572.CEL.gz 
        10.38333         10.31405         10.25655

qq显示，这个数据比较接近于正态分布，但是t检验发现，这个基因在两组之间并无差异，如下所示：

t.test(e[g==1],e[g==0])$p.value

结果如下所示：

> t.test(e[g==1],e[g==0])$p.value
[1] 0.779303

上面是只计算1个基因的情况，如果要计算所有的基因，那么就需要对每个基因进行t检验，我们来构建一个函数，如下所示：

myttest <- function(x) t.test(x[g==1],x[g==0],var.equal=TRUE)$p.value
pvals <- apply(geneExpression,1,myttest)
head(pvals)

结果如下所示：

> head(pvals)
 1007_s_at    1053_at     117_at     121_at  1255_g_at    1294_at 
0.04553344 0.03370683 0.13604026 0.59413846 0.96849102 0.08489586

现在我们来看一下，p值小于0.05的基因有多少个，如下所示：

sum(pvals<0.05)

结果如下所示：

> sum(pvals<0.05)
[1] 1383

也就是说，有1383个基因有差异。

但是，正如我们随后即将详细描述的那样，对于这个结果，我们要谨慎对待。因为这个结果是我们经过8000多次比较后得到的结果（数据集里一共有个8793基因）。如果我们对一批随机数据执行相同的计算过程，在零假设成立的前提下，我们也能计算出一批p值，并且这批p值小于0.05的数目为：

set.seed(1)
m <- nrow(geneExpression)
n <- ncol(geneExpression)
randomData <- matrix(rnorm(n*m),m,n)
nullpvals <- apply(randomData,1,myttest)
sum(nullpvals<0.05)

结果如下所示：

> sum(nullpvals<0.05)
[1] 419

随后，我们会解释一下这个数值是怎么一回事，实际上，这个数字大概就是8192*0.05的结果（大约为409.6），后面会详细介绍这个原理。

加快t检验的计算

在我们开始之前，我们需要知道，前面我们构建了一个做t检验的函数，虽然它能计算出相应的p值，但是这种方法比较低效。我们还有其它更快的方法实现这一目的，在Bioconductor项目中就有这样的函数，此时我们要用到rafalib包中的install_bioc()函数，如下所示：

install_bioc("genefilter")

现在我们使用rowttests()函数来计算一下p值，如下所示：

library(genefilter)
results <- rowttests(geneExpression,factor(g))
max(abs(pvals-results$p))
sum(results$p<0.05)

结果如下所示：

> max(abs(pvals-results$p))
[1] 6.528111e-14
> sum(results$p<0.05)
[1] 1383

p值计算练习

为了更好地理解p值是一个随机变量，现在我们再来看一些案例，如下所示：

set.seed(1)
library(downloader)
url = "https://raw.githubusercontent.com/genomicsclass/dagdata/master/inst/extda\
ta/femaleControlsPopulation.csv"
filename = "femaleControlsPopulation.csv"
if (!file.exists(filename)) download(url,destfile=filename)
population = read.csv(filename)

pvals <- replicate(1000,{
  control = sample(population,12)
  treatment = sample(population,12)
  t.test(treatment,control)$p.val})
head(pvals)
hist(pvals)

结果如下所示：

> head(pvals)
[1] 0.76454370 0.60237592 0.49679830 0.23518016 0.20882028 0.06812489

现在回答以下几个问题：

问题一：小于0.05的p值所占的比例是多少？

答：4.1%，代码如下所示：

> sum(pvals < 0.05)/length(pvals)*100
[1] 4.1

问题二：小于0.01的p值所占的比例是多少？

答：0.8%，如下所示：

> sum(pvals < 0.01)/length(pvals)*100
[1] 0.8

问题三：假设你正在测试20种饮食对小鼠体重影响的效应。对于20种饮食中的每一种，我们都使用10只对照小鼠和10只实验小鼠。在零假设成立时，也就是说饮食对小鼠的体重没有影响，另外，小鼠的体重是服从正态分布的（均值为30g，SD为2g）。我们对其中的一项进行Monte Carlo模拟：

cases = rnorm(10,30,2)
controls = rnorm(10,30,2)
t.test(cases,controls)

结果如下所示：

> cases = rnorm(10,30,2)
> controls = rnorm(10,30,2)
> t.test(cases,controls)

	Welch Two Sample t-test

data:  cases and controls
t = 0.16473, df = 17.934, p-value = 0.871
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.708669  1.999327
sample estimates:
mean of x mean of y 
 30.23172  30.08639

问题：现在进行Monte Carlo模拟，模拟所有20种包含的实验结果，如果我们将种子设置为set.seed(100)，那么有多少值低于0.05？

答案是：5%，如下所示：

set.seed(100)
pvals <- replicate(20,{
  cases = rnorm(10,30,2)
  controls = rnorm(10,30,2)
  t.test(cases,controls)$p.val})
sum(pvals<0.05)/length(pvals)*100

结果如下：

> sum(pvals<0.05)/length(pvals)*100
[1] 5

问题四：现在我们通过模拟数据了解了p值的分布情况，也就是说了解了小于0.05的p值的比例。在问题三中，我们一次运行了20种饮食的实验。现在我们运行1000次实验，并且每次保存下那些p值小于0.05的数目。

操作如下：

随机种子数设为100，即set.seed(100)，运行问题三中的代码1000次，并保存p值小于0.05的次数，计算一下这些数字的平均值，这就是我们拒绝零假设为真时的次数，如下所示：

set.seed(100)
pvals <- replicate(1000,{
  cases = rnorm(10,30,2)
  controls = rnorm(10,30,2)
  t.test(cases,controls)$p.val})
sum(pvals<0.05)
sum(pvals<0.05)/length(pvals)*100

结果如下所示：

> sum(pvals<0.05)
[1] 40
> sum(pvals<0.05)/length(pvals)*100
[1] 4

问题五：这说明，就平均数而言，即使所有的饮食对小鼠的体重都没有影响，我们还能得到一些小于0.05的p值。使用相同的模拟数据来做1000次计算，我们能得到多少假阳性的比例？

答：大概是5%。